#SLB Nest SDM 

sapply(c("raster","rgdal","sp","sf","ggplot2", "tidyverse", "fasterize"), library, character.only=T)


#Get files ready ####
#setwd("") #insert working directory with files

#Buffer area (Blue Mountains extent)
BM_buff <- readOGR("BM_3km_buff2.shp")
crs(BM_buff)
proj4string(BM_buff) <-CRS('+proj=merc +zone=56 +south +datum=WGS84 +units=m +no_defs')
plot(BM_buff, add=TRUE)

#Vegetation
veg <- raster("VegMap4.tif")
veg
plot(veg)
proj4string(veg) <-CRS('+proj=merc +zone=56 +south +datum=WGS84 +units=m +no_defs')
crs(veg)

#Key for vegetation types
vegdata <- read.csv("VegKey_VM4.csv")
View(vegdata)

#Reclassify vegetation layer to four vegetation types
rc <- as.matrix(data.frame(from=vegdata$Val, to=vegdata$NewNo))
new_veg <- reclassify(veg, rc)
new_veg
plot(new_veg)
plot(BM_buff, add=TRUE)

#Slope
#Export from ArcGIS in same coordinate system as veg map (WGS mercator)
slope <- raster("Slope_DEM_10m2_WGS_merc.tif")
slope
plot(slope)
crs(slope)
plot(BM_buff, add=TRUE)

#Fire
fire <- raster("FESM_1.tif")
crs(fire)
plot(fire)
plot(BM_buff, add=TRUE)


#Split vegetation layer for ArcGIS analysis #####
new_veg2 <- mask(new_veg, BM_buff)
plot(new_veg2)
plot(BM_buff, add=TRUE)
rat <- data.frame(
  ID = 1:4,
  vegType = c("Drysclerophyll", "Rainforest", "Wetsclerophyll", "Other")
)
levels(new_veg2) <- rat
veg4 <- for(i in 1:4){
  assign(rat$vegType[i],mask(new_veg2,new_veg2 != i, maskvalue=1))
}
plot(Rainforest)
plot(Drysclerophyll)
plot(BM_buff, add=TRUE)
plot(Wetsclerophyll)
plot(Other)
proj4string(Rainforest) <-CRS('+proj=merc +zone=56 +south +datum=WGS84 +units=m +no_defs')
proj4string(Drysclerophyll) <-CRS('+proj=merc +zone=56 +south +datum=WGS84 +units=m +no_defs')
proj4string(Wetsclerophyll) <-CRS('+proj=merc +zone=56 +south +datum=WGS84 +units=m +no_defs')
proj4string(Other) <-CRS('+proj=merc +zone=56 +south +datum=WGS84 +units=m +no_defs')
writeRaster(Rainforest, "Rainforest.tif")
writeRaster(Drysclerophyll, "Drysclerophyll.tif")
writeRaster(Wetsclerophyll, "Wetsclerophyll.tif")
writeRaster(Other, "Otherveg.tif")


#Split fire raster for ArcGIS analysis #####
fire2 <- mask(fire, BM_buff)
plot(fire2)
#check number of levels
barplot(fire2) #only 5, fesm is missing
fire2
#Check which levels are present in map
plot(fire2, col = c("white", "green", "yellow", "orange", "red", "black"))
fire2
#remove fesm - level 1
fire.rc <- data.frame(
  Val = c(0, 1, 2, 3, 4, 5),
  newVal = c(1, 1, 2, 3, 4, 5)
)
View(fire.rc)
rc.f <- as.matrix(data.frame(from=fire.rc$Val, to=fire.rc$newVal))
rc.f
new_fire <- reclassify(fire2, rc.f)
new_fire
plot(new_fire, col = c("white", "yellow", "orange", "red", "black"))

new_fire2 <- mask(new_fire, BM_buff)
plot(new_fire2, col = c("white", "gold", "orange", "red", "black"))
writeRaster(new_fire2, "Fire_noFESM.tif")

rat2 <- data.frame(
  ID = 1:5,
  fire_sev = c("Unburnt","Low","Moderate","High","Extreme")
)
View(rat2)
levels(new_fire2) <- rat2
new_fire2
fesm <- for(i in 1:5){
  assign(rat2$fire_sev[i],raster::mask(new_fire2,new_fire2 != i, maskvalue=1))
}
plot(Unburnt)
plot(Low)
plot(Moderate)
plot(High)
plot(Extreme)
writeRaster(Unburnt, "Unburnt.tif")
writeRaster(Low, "LowFire.tif")
writeRaster(Moderate, "ModFire.tif")
writeRaster(High, "HighFire.tif")
writeRaster(Extreme, "Extreme.tif")


#Prepare data for SDM #####

#Vegetation
veg2 <- resample(new_veg, slope, method = 'ngb')

#Restrict layers to Blue Mountains only
veg3 <- raster::mask(veg2, BM_buff)
plot(veg3)
slope3 <- mask(slope, BM_buff)
plot(slope3)

#Add nest locations
nest_locs <- read.csv("Nest_data.csv")
View(nest_locs)

#Prepare nest points to get pres and abs points
pres_locs <- nest_locs[nest_locs$PA==1,]
View(pres_locs)
abs_locs <- nest_locs[nest_locs$PA==0,]
View(abs_locs)


#### Run SDM ####
library(glm2)
library(dismo)
library(maptools)
library(sp)

#plot nest points for visualising
data(wrld_simpl)
plot(wrld_simpl, xlim=c(148, 152), ylim=c(-35, -32), axes=TRUE, col="light yellow")
box()
points(nest_locs$Longitude, nest_locs$Latitude, col="red", pch=20, cex=0.75)

coordinates(pres_locs) <- ~Longitude+Latitude
crs(pres_locs) <- crs(wrld_simpl)
class(pres_locs)

coordinates(abs_locs) <- ~Longitude+Latitude
crs(abs_locs) <- crs(wrld_simpl)
class(abs_locs)

crs(pres_locs)
crs(abs_locs)

#transform to match other layers
pres_locs <- spTransform(pres_locs, crs(slope))
abs_locs <- spTransform(abs_locs, crs(slope))
plot(slope3)
plot(veg3)
plot(abs_locs, col = "red", pch=20, cex=0.75, add=T)
plot(pres_locs, col = "blue", pch=20, cex=0.75, add=T)

#Put data together for SDM
predictors <- stack(veg3, slope3)
presvals <- raster::extract(predictors, pres_locs)
absvals <- raster::extract(predictors, abs_locs)

#create data frame for GLM
train <- rbind(pres_locs, abs_locs)
pb_train <- c(rep(1, nrow(pres_locs)), rep(0, nrow(abs_locs))) #creates string of presence/absence to match data
envtrain <- raster::extract(predictors, train)
envtrain <- data.frame(cbind(pa=pb_train, envtrain)) #creates data fram with nest presence/absence and spatial variables
View(envtrain)
envtrain[,'VegMap4'] = factor(envtrain[,'VegMap4'], levels=1:4)

#Run GLM
gm2 <- glm(pa ~ VegMap4 + Slope_DEM_10m2_WGS_merc, family = binomial(link="logit"), data=envtrain)
summary(gm2)
#predict suitability based on GLM
pg2 <- predict(predictors, gm2)
pg2
writeRaster(pg2, "suitability.tif")
plot(pg2, main = 'Overall suitability')

#Create categories of suitability using natural breaks
library(BAMMtools)
pg2_vals <- rasterToPoints(pg2)
head(pg2_vals)
class(pg2_vals)
pg_dat <- as.data.frame(pg2_vals)
head(pg_dat)
View(pg_dat)
sdm_vals <- pg_dat$layer #check that this variable is named correctly - it should be the suitability
head(sdm_vals)

#Break into six categories
getJenksBreaks(sdm_vals, 5, subset = 10000)
#-18.3435040 -13.1810379  -0.6241069   0.5543447   3.9167016

#find min and max values of sdm
pg2
#-18.344, 6.042365

#reclassify into four categories
reclass <- c(-19, -13.1810379  , 1,
             -13.1810379  , -0.6241069, 2,
             -0.6241069, 0.5543447, 3,
             0.5543447, 7, 4)
reclass_m <- matrix(reclass, ncol=3, byrow=TRUE)
reclass_m
nest_sdm_class6.4 <- reclassify(pg2, reclass_m) #takes 2 minutes
barplot(nest_sdm_class6.4)
plot(nest_sdm_class6.4)
nest_sdm_class6.4
writeRaster(nest_sdm_class6.4, "sdm_class.tif")

#split raster into categories of suitability
#split raster
rat <- data.frame(
  ID = 1:4,
  Suitability = c("Unsuitable", "LowSuit", "ModSuit", "HighSuit")
)
View(rat)
levels(nest_sdm_class6.4) <- rat
nest_sdm_class6.4
sdm_class <- for(i in 1:5){
  assign(rat$Suitability[i],mask(nest_sdm_class6.4,nest_sdm_class6.4 != i, maskvalue=1))
}
plot(Unsuitable)
plot(LowSuit)
plot(ModSuit)
plot(HighSuit)
writeRaster(Unsuitable, "Unsuitable.tif")
writeRaster(LowSuit, "LowSuit.tif")
writeRaster(ModSuit, "ModSuit.tif")
writeRaster(HighSuit, "HighSuit.tif")